Flagellum-driven cargoes: Influence of cargo size and the flagellum-cargo attachment geometry

The beating of cilia and flagella, which relies on an efficient conversion of energy from ATP-hydrolysis into mechanical work, offers a promising way to propel synthetic cargoes. Recent experimental realizations of such micro-swimmers, in which micron-sized beads are propelled by isolated and demembranated flagella from the green algae Chlamydomonas reinhardtii (C. reinhardtii), revealed a variety of propulsion modes, depending in particular on the calcium concentration. Here, we investigate theoretically and numerically the propulsion of a bead as a function of the flagellar waveform and the attachment geometries between the bead and the flagellum. To this end, we take advantage of the low Reynolds number of the fluid flows generated by the micro-swimmer, which allows us to neglect fluid inertia. By describing the flagellar waveform as a superposition of a static component and a propagating wave, and using resistive-force theory, we show that the asymmetric sideways attachment of the flagellum to the bead makes a contribution to the rotational velocity of the micro-swimmer that is comparable to the contribution caused by the static component of the flagellar waveform. Remarkably, our analysis reveals the existence of a counter-intuitive propulsion regime in which an increase in the size of the cargo, and hence its drag, leads to an increase in some components of the velocity of the bead. Finally, we discuss the relevance of the uncovered mechanisms for the fabrication of synthetic, bio-actuated medical micro-robots for targeted drug delivery.


Introduction
Cilia and flagella are fundamental units of motion in various biological systems. These hairlike cellular protrusions share a common conserved 9+2 microtubule-based structure, and beat to accomplish a variety of biological tasks. Examples are surface fluid flows generated by ciliary carpet in the human respiratory tract to remove pollutants [1], cilia-driven cerebrospinal fluid transport in the mammalian brain to deliver nutrients and important signaling molecules [2], ciliary flow in the Fallopian tube to assist sperm transport to the fertilization site [3], and propulsion of green algae C. reinhardtii that swims by breaststroke-like motion of its two flagella [4][5][6].
In recent years, there has been a great interest in the field of targeted drug delivery and assisted fertilization to integrate cilia and flagella as efficient energy conversion modules into bio-compatible micro-swimmers. Autonomous flagella-driven motility of various biological species, mainly E. coli and sperm, are utilized as bio-actuators to provide an efficient cargo transport [7][8][9][10][11][12][13]. More recently, in the experiments by Ahmad et al. [14], axonemally-driven cargoes are fabricated by integration of isolated and demembranated flagella from C. reinhardtii (known as axonemes) with micron-sized beads (see Fig 1 and S1-S2 Videos). These ATPreactivated axonemes, with a length of approximately 10 μm, beat with an ATP-dependent frequency [15][16][17] and have an asymmetric waveform that can best be described as a base-to-tip traveling wave component superimposed on a circular arc with mean curvature of about -0.2 μm −1 . The static component of the axonemal curvature leads to a curved swimming trajectory of the micro-swimmer (see Fig 1D) [18][19][20][21]. In comparison, when the static curvature is strongly reduced, a micro-swimmer swims along an essentially straight path [14,[22][23][24]. Importantly, the static curvature of axonemes is highly dependent on the calcium concentration. Namely, increasing the calcium concentration beyond 0.05 mM reduces the static curvature of axonemes by one order of magnitude [14,23,24], thereby inducing a transition from circular to straight swimming trajectories of axonemally-propelled beads. In addition to the flagellar waveform, the axoneme-bead attachment geometry also plays a critical role in the cargo propulsion speed. As emphasized in Ref. [14], axonemes can attach to the bead symmetrically, with their tangent vector at the contact point passing through the bead center, or asymmetrically. Due to the limitations of 2D microscopy in Ref. [14], it was not experimentally possible to distinguish between these two types of bead-axoneme attachment and 3D microscopy techniques [25] are required to quantify the effect of attachment geometry on the propulsion speed. This symmetric versus asymmetric attachment has consequences on cargo propulsion dynamics and investigating this effect theoretically and numerically is the main focus of the present work.
Here, we investigate the effect of (i) various flagellar wave components, (ii) the size of the cargo (the bead), and (iii) the symmetric versus asymmetric flagellum-bead attachment on the swimming dynamics of a bead that is propelled by a model flagellum. We restrict ourselves to two-dimensional motion, which captures most of the experimental results in Ref. [14]. We use an approximate description of the flagellum waveform as a combination of a static curvature and a traveling wave component, and use resistive-force theory (RFT) [26,27] to obtain analytical expressions for the translational and rotational velocities of a flagellum-propelled bead in the limit of small amplitude of curvature waves. We compare the resulting expressions to the results of simulations of the swimming trajectories. Our analysis reveals a surprising nonmonotonic behavior of the mean translational and rotational velocities of the axonemallydriven bead as a function of the bead radius. Finally, our analysis shows that for a freely swimming axoneme, which rotates predominantly with its static component of the axonemal waveform, sideways bead attachment is sufficient to generate mean rotational velocities comparable to the rotation rates induced by the static curvature. This paper is structured as follows: first we briefly describe RFT, which we use to calculate the propulsion speed of the micro-swimmer as a function of the cargo size. In this approach, the specific details of the bead-axoneme attachment geometry are taken into account in the drag matrices of the bead and the axoneme.
Next, we present our analytical approximations and numerical simulations to show the effect of the cargo size and of the symmetric versus asymmetric bead-axoneme attachment.

RFT and calculations of mean translational and rotational velocities
The fluid flow generated by the swimming of small objects is characterized by very small Reynolds numbers. In this regime, viscous forces dominate over inertia and non-reciprocal motion is necessary to break the time-symmetry and generate propulsion (scallop theorem) [28,29]. The micro-swimmer in our system consists of an axoneme (a filament of characteristic length L * 10 μm and radius 0.1 μm), which is attached at one end to a micron-sized bead and swims in an aqueous solution of viscosity μ = 10 −3 Pa s and density ρ = 10 3 kg m −3 . Given the characteristic axonemal wave velocity V = λ/T � 0.5 mm s −1 (calculated for a typical from green algae C. reinhardtii is attached to a 1 micron-sized bead. The axoneme is reactivated with 1 mM ATP and beats at around 110 Hz (see S1 Video). B) Over time, the position of the center of the bead, its evolution represented by the blue curve, is propelled on a helical-like trajectory (see S2 Video). (C) The definition of the laboratory and the swimmer-fixed frame. As the flagellum beats, the micro-swimmer swims counterclockwise (CCW) in the microscope's field of view effectively in 2D. D) The traces of the basal (yellow line) and distal tip (cyan line) of the flagellum tracked for 198 sec. (E) Curvature waves initiate at the basal end of axoneme (s = 0) which is attached to the bead, and propagate toward the distal tip (s = L) at the frequency of about 110 Hz. F) The axonemal shapes averaged over one beat cycle results in a circular arc with mean curvature of about -0.2 μm −1 . This static component of the axonemal curvature results in a curved swimming trajectory and in the absence of this component, the bead is propelled on a straight trajectory. The experimental techniques used to record the motion of the beads and the flagella are described in Ref. [14]. Please note that in this experiment, the bead-axoneme attachment appears to be asymmetric but 3D microscopy techniques are required to distinguish a symmetric versus an asymmetric attachment.
https://doi.org/10.1371/journal.pone.0279940.g001 axonemal beat frequency of 50 Hz and a wavelength which is comparable to the axonemal contour length L), the Reynolds number Re = ρLV/μ is small, no larger than * 0.005. In this physical regime, Newton's laws then consist of an instantaneous balance between external and fluid forces and torques exerted on the swimmer, i.e. F ext + F fluid = 0 and τ ext + τ fluid = 0. The force F fluid and torque τ fluid exerted by the fluid on the axoneme-bead swimmer can be written as: where F bead and τ Bead are the hydrodynamic drag force and torque acting on the bead, and the integrals over the contour length L of the axoneme calculate the total hydrodynamic force and torque exerted by the fluid on the axoneme. The bead is propelled by the oscillatory shape deformations of the ATP-reactivated axoneme. At any given time, we consider axoneme-bead swimmer as a solid body with translational and rotational velocities U(t) and O(t) to be determined as explained below. F fluid and τ fluid can be separated into propulsive part due to the relative shape deformations of the axoneme in the body-fixed frame and the drag part [30]: where the 6 × 6 geometry-dependent drag matrix D is symmetric and non-singular (invertible) and is composed of drag matrix of the axoneme D A and drag matrix of the bead D B . For a freely swimming axoneme-bead, which experiences no external forces and torques, F fluid and τ fluid must vanish. As explained in the introduction, we restrict ourselves to 2-dimensional motion, which describes most of the experimental work described by Ref. [14]. In 2-dimensions, D is reduced to a 3 × 3 matrix and Eq 3 can be reformulated as: which we use to calculate translational and rotational velocities of the swimmer after determining the drag matrices D A and D B , and the propulsive forces and torque ðF prop in the body-fixed frame by selecting the basal end of the axoneme (bead-axoneme contact point) as the origin of the swimmer-fixed frame. As shown in Figs 1C and 3A, we define the local tangent vector at contour length s = 0 as X-direction, the local normal vector n as the Y-direction, and assume that z and Z are parallel. Here (x, y, z) denote an orthogonal lab-frame basis. We define θ 0 (t) = θ(s = 0, t) as the angle between x and X which gives the velocity of the bead in the laboratory frame as U BeadÀ Lab for a given beating pattern of axoneme in the body-fixed frame, we used the classical framework of resistive-force theory (RFT), which neglects longrange hydrodynamic interactions between different parts of the flagellum as well as the inter-flagella interactions [26,27]. In this theory, the flagellum is discretized as a set of small rod-like segments moving with velocity u 0 (s, t) in the body-frame, as illustrated in Fig 2. The propulsive force F prop is proportional to the local center-line velocity components of each segment in parallel and perpendicular directions: where u 0 k and u 0 ? are the projections of the local velocity on the directions parallel and perpendicular to the axoneme. The friction coefficients in parallel and perpendicular directions are defined as z k = 4πμ/(ln(2L/a) + 0.5) and z ? = 2z k [26], respectively. This anisotropic hydrodynamic friction experienced by a cylindrical segment is the basis of propulsion by a flagellum. For a water-like environment with dynamic viscosity μ = 10 −3 Pa s, we obtain z k * 2.1 pN msec/μm 2 . As z k 6 ¼ z ? , Eq (5) implies that the resulting velocity is not parallel to the propulsive force F prop . In the following, we introduce the two dimensionless quantities: The value of η will be fixed to 1/2 in the rest of the text [26]. The value of z 0 ? is determined by the geometry of the axoneme. We take for the axoneme radius a realistic value of a = 0.1 μm and a contour length of L = 10 μm.
Here is a brief summary of the steps in the RFT analysis: First, we translate and rotate the axoneme such that the basal end is at position (0, 0) and the local tangent vector at s = 0 at any time t is along the x-axis (see Fig 1C). In this way, we lose the orientation information of the axoneme at all the time points except for the initial configuration at time t = 0. Second, we calculate propulsive forces and torque in the body-frame using RFT (Eq 5), and then use Eq 4 to obtain translational velocities U x , U y as well as rotational velocity O z of the axoneme. Now the infinitesimal rotational matrix can be expressed as which we use to update the rotation matrix as Γ(t + dt) = Γ(t)dΓ(t), considering Γ(t = 0) to be the unity matrix. Having the rotation matrix at time t, we obtain the configuration of the axoneme at time t from its shape at the body-fixed frame by multiplying the rotation matrix as r lab-frame (s, t) = Γ(t)r body-frame (s, t), which can then be compared with the experimental data. Please note that r lab-frame (s, t) = (X lab-frame (s, t), Y lab-frame (s, t), 1) is an input from experimental data presenting the beating patterns in the laboratory frame.
S2 Fig in S1 File shows a comparison between the rotational and translational velocities measured directly with the experimental data presented in Fig 1 and the results obtained with RFT using the experimental beat pattern as input (as explained above). This comparison shows a semi-quantitative agreement, therefore justifying our analysis in the framework of RFT in section 3.2.

Drag matrix of a bead in 2D
Let us consider the two-dimensional geometry defined in Fig 3A. Note that the origin of the swimmer-fixed frame is not at the bead center; rather it is selected to be at the bead-axoneme contact point. In general, the tangent vector at position s = 0 of the axoneme, which defines the X-axis, does not pass through the bead center located at ( . This asymmetric bead-axoneme attachment is also observed in the experiments, as shown in Fig 1A. We emphasize that the drag force is actually a distributed force, given by d f = σ.d A, applied at the surface of the sphere, but symmetry implies that drag force effectively acts on the bead center. We define the translational and rotational friction coefficients of the bead as ν T = 6πα t μR and ν R = 8πα r μR 3 , where μ = 10 −3 Pa s is the dynamic viscosity of water and factors α t = 1/(1 − 9(R/h)/16 + (R/h) 3 /8) and α r = 1/(1−(R/h) 3 /8) are corrections due to the fact that axonemal-based bead propulsion occurs in the vicinity of a substrate [31]. Here R is the bead radius (*0.5 μm), and h is the distance of the center of the bead to the substrate. Assuming R/h * 1, we obtain α t = 16/9 and α r = 8/7. We now look at each component of velocity and ask what force do we need to apply to counteract the viscous force and torque?
i. Translation in X-direction. In this case, we have (U x , U y , O z ) = (1, 0, 0) as shown in Fig 3B. We need to apply a force in +X-direction to counteract the drag force as: But we must also apply a torque in +Z-direction for the case illustrated in Fig 3B (where Y B < 0) to prevent rotation from occurring: so we get: ii. Translation in Y-direction, see Fig 3C, corresponding to (U x , U y , O z ) = (0, 1, 0). We have F x = 0 and F y = + 6πα t μR = ν T . Note that we need to apply a negative torque, and since X B < 0, we have τ z = + X B ν T which gives: iii. Rotation around Z-direction, see Fig 3D, corresponding to (U x , U y , O z ) = (0, 0, 1). Before looking at the forces, let us examine the motion. The rotation O 0 z ¼ 1 of the bead center around the origin O also generates translational velocity x > 0 and U 0 y < 0 which is consistent. Around the center of the bead, drag exerts force and torque F D and τ D , as depicted in Fig 3D: To counteract the drag force, we must apply: Now we combine parts (i), (ii) and (iii) to obtain: For the special case of a symmetric attachment, with the center of the bead at (X B , Y B ) = (−R, 0), Eq 14 simplifies to: Note that Eqs 14 and 15 present the forces and torque exerted by the bead on the fluid which has opposite sign of the forces generated by the fluid on the bead, so the drag matrix of the bead D B is given by: The general form of the drag matrix in 3D is derived in the supplemental information.

Drag matrix of an axoneme in 2D
Since the axoneme beats over time, its drag matrix, which relates force and velocity, is timedependent: Assume the motion is 2D, and consider the axonemal shapes at successive time points are given in the body frame as r body-frame (s, t) = (X body-frame (s, t), Y body-frame (s, t), 1). This can be either an input from the experimental data (see Fig 1F) or a predefined waveform of a model axoneme, as introduced in Section 3.1. The experimental shapes shown in Fig 1 were recorded with a high time resolution of 1000 Hz [14], and translated and rotated such that the tangent vector at s = 0 is along the X-axis. Factoring out an overall rotation and translation in the laboratory frame allows us to focus on the shape deformation of the axoneme.
To obtain the elements of the drag matrix for a given axonemal shape at time t, we work in the framework of RFT and follow the same procedure as in the case of a bead, as described in Sec. 2.2: i. Global translation of the axoneme in X-direction. In this case, we have (U x , U y , O z ) = (1, 0, 0), which using Eq 5 and the given axonemal shape at time t, we first obtain the tangential and perpendicular velocity components for each cylindrical segment of the axoneme. Second, in the framework of RFT, we calculate the corresponding elemental force dF = (dF X (s, t), dF Y (s, t)) and torque dτ Z (s, t) = r body-frame (s, t) × dF exerted by each cylindrical segment of the axoneme on the fluid to counteract the drag force. Third, to obtain the drag elements a 11 , a 21 and a 31 , we integrate the elemental force and torque over the whole contour length of axoneme to calculate the total force and torque exerted by the axoneme on the fluid to counteract the drag. The fluid drag force exerted on the axoneme has an opposite sign, thus: ii. Global translation of the axoneme in Y-direction. In this case, we have (U x , U y , O z ) = (0, 1, 0). The matrix elements a 12 , a 22 and a 32 are obtained as in (i).
iii. Global rotation of the axoneme around Z-direction, corresponding to (U x , U y , O z ) = (0, 0, 1). We note that the rotation of the axoneme around the origin O (see Fig 3A) with O = (0, 0, O z ), also generates translational velocity components as O × r body-frame (s, t) = (−Y body-frame (s, t), X body-frame (s, t)), that should be taken into account while using Eq 5 to obtain the force and the torque that the axoneme exerts on the fluid to counteract the drag. Similar to Eq 18, we can now calculate the drag elements a 13 , a 23 and a 33 .

Analytical approximations of rotational and translational velocities of an axonemally-propelled bead
The waveform of the axoneme is complex, and involves a combination of several components [18,20,22,32]. In this first subsection, we begin by discussing a simplified waveform, in order to understand the elementary aspects of the propulsion of the bead. In practice, we approximate the waveform of the axoneme as a superposition of traveling wave component, with amplitude C 1 , and a circular arc with mean curvature C 0 [18,22]: where ω 0 = 2πf 0 , k = 2π/λ is the wave number. For our exemplary axoneme in Fig 1A, following the method described in Ref. [16], we calculate the wavelength to be λ * 11.34 μm, which is *34% larger than the axonemal contour length L * 8 μm. The approximate waveform given by Eq 19 allows us to obtain explicit expressions for the propulsion velocity of the cargo. This expression, however, neglects a small backwards wave component, propagating from tipto-base of the form C 0 1 cosðo 0 t þ ksÞ [22], along with components with wave numbers equal to n × k, where n is an integer > 1. The results of our analysis of beating axonemes [22] show that back-propagating wave component is about 5-10 times smaller than the main base-to-tip wave. For simplicity, this small component is neglected in this subsection. The analysis presented in Subsection 3.2, based on numerical simulations with the precise waveform of the axonemes determined experimentally [14], qualitatively validates the approach presented here.
To obtain analytical expressions for the mean translational and rotational velocities of the swimmer, we used RFT, as presented in the Materials and Methods section. As a further simplification, we neglect in this section the difference between the wavelength of the beat pattern, λ, and the size of the axoneme, L, and assume in the following L = λ. We determine the propulsion velocity up to the first order in C 0 , and the second order in C 1 .

Symmetric bead-axoneme attachment.
Let us first consider the example of an axoneme attached symmetrically from the basal side to a bead, so that the tangent vector at s = 0 passes through the bead center (X B = −R, Y B = 0, see Fig 3B). We determine the dependence of the translational and rotational velocities of the swimmer on the dimensionless bead radius r = R/L, with the simplified waveform of the axoneme given by Eq 19 and impose the force-free and torque-free conditions in 2D. The drag matrix of the bead is given by Eq 16, with Y B = 0 and X B = −R. The drag matrix of the axoneme is calculated as described in Sec. 2.2.
We approximate analytically the averaged angular and linear velocities in the swimmerfixed frame, as defined in Fig 1C, and we determine the propulsion velocities up to first order in C 0 , and to second order in C 1 : The functions β 1 , β 2 and β 3 depend on the dimensionless quantities z 0 ? and η, defined by Eq 6, and on r = R/L. The explicit expressions are presented in the supplemental information, Eqs. S16-S18. The results in the absence of the bead simplifies to: which are previously discussed in Refs. [19,20]. The corresponding dependence on C 0 and C 1 is shown in Fig 4A-4F as full lines. We note that Eqs 20-22, in the absence of intrinsic curvature C 0 = 0, predict that the axoneme swims in a straight path with hU y i = 0, hU x i proportional to the square of traveling wave component C 1 [33][34][35] and the mean rotational velocity hO z i vanishes (see the solid red line in Fig 4A and S3A Fig in S1 File).
In order to verify the quality of our analytical approximations, we also determined the motion of the swimmer numerically using RFT, starting from the simplified waveform given in Eq 19 and r = 0. The corresponding results are shown by the circular symbols in Fig 4A-4F. The comparison between numerical simulations and the full analytical approximations presented in Eqs. S16-S18, shows a very good agreement at small values of C 0 and C 1 , with deviations at larger values. In addition, three exemplary trajectories (r = 0.1), determined from RFT, are shown in Fig 4G-4I. The corresponding averaged rotational velocity hO z i of the model swimmer is proportional to the square of the traveling wave component C 1 .
To investigate the dependency of the mean translational and rotational velocities of our model swimmer on the bead size, r, Fig 5 shows the dependence of hO z i/ω 0 (A and D), hU x i/ (Lω 0 ) (B and E) and hU y i/(Lω 0 ) (C and F), predicted by Eqs 20-22; see the full lines. We also performed numerical simulations at different values of the bead radii, see the circular symbols. As shown in Fig 5, there is a very good agreement between our numerical simulations and analytical approximations at small values of C 1 (panels A-C) but deviations appear at larger values (panels D-F). Remarkably, while hU x i decreases monotonically with the bead radius r = R/L, both hO z i and hU y i exhibit a non-monotonic dependence. This behavior is counter-intuitive: the drag exerted by the fluid on the sphere increases with the size of the bead, which in turn increases dissipation. Based on this general remark, one expects the velocity of the swimmer to go down when r increases. We nevertheless notice that not all three components of velocity may increase when r increases.
To gain more insight on this anomalous behavior, we determined the asymptotic expressions of Eqs 20-22 in two opposite limits of small and large bead radii. The corresponding dependence is shown by the dashed lines in Fig 5. In the limit of very large bead radius, R � L (small 1/r), we obtain a dependence of hU x i and hU y i as r −1 , and of hO z i as r −2 (up to the higher order corrections): þ720ð3Z À 1Þr À 1 À 96p 2 ðZð9r À 1 þ 30Þ À 11r À 1 À 30ÞÞ: In the opposite limit of small r (i.e. R � L), up to the second order in r, we obtain with the realistic value of z 0 ? � 4:33: hO z i o 0 � À 0:42C 0 C 2 1 ð1 À 19:15r þ 513:77r 2 Þ; ð29Þ

PLOS ONE
Flagellum-driven cargoes The black dashed lines (with and without stars in pink color) in Fig 5 show the corresponding asymptotic behavior. The limiting behavior is valid only for very small values of r. It nevertheless qualitatively captures the non-monotonous trend, observed at intermediate values of r. The transition from the r 2 dependence at small values of bead radius to the r −2 -trend at large values of r provides a qualitative explanation for the non-monotonous behavior, observed for hO z i and hU y i.

The sideways bead-axoneme attachment contributes to the rotational velocity of the swimmer.
In the experiments in Ref. [14], it was frequently observed that the bead-axoneme attachment was asymmetric, i.e. the tangent vector of the axoneme at s = 0 does not pass through the bead center. This case is schematically illustrated in Fig 3A, where it results in a value of Y B 6 ¼ 0. Interestingly, our analytical approximations and simulations show that this asymmetric bead-axoneme attachment is enough to rotate the axoneme, so the presence of the static curvature or the second harmonic is not necessary for rotation to occur.
For this analysis, we consider the 2D geometry where the center of the bead is at position X B and Y B , measured with respect to the coordinate system defined at the bead-axoneme contact point (Figs 1C and 3A). We will use here the dimensionless coordinates . The drag matrix of the bead is given by Eq 16, where we specify the value of Y B 6 ¼ 0 corresponding to the asymmetric attachment. The beating of the axoneme is described by Eq 19 in terms of a traveling wave component C 1 and intrinsic curvature C 0 . Similar to the case of a symmetric bead-axoneme attachment in Section 3.1.1, we calculate the mean rotational and translational velocities of an axonemally-driven bead by combining the drag matrix of the bead and the axoneme (see Materials and methods).
The results up to the leading order in C 0 and C 1 can be expressed as: where η and z 0 ? are defined by Eq 6, r = R/L, and y b = Y b /L. When the attachment is symmetric (y B = 0), symmetry considerations impose that: which expresses the fact that in the absence of C 0 , both hO z i and hU y i are zero. This is also consistent with previous expressions of hO z i and hU y i when the attachment is symmetric, see Eqs 26,28,29 and 31. In contrast, when the attachment is not symmetric (y B 6 ¼ 0), the coefficients a 1 ðZ; z 0 ? ; y B ; rÞ and a 3 ðZ; z 0 ? ; y B ; rÞ become non zero, which implies that the system rotates (hO z i6 ¼0) and has a nonzero velocity component hU y i6 ¼0 even when C 0 = 0. This clearly shows the importance of the asymmetric bead-axoneme attachment. This is illustrated by Fig 6, which shows the translational and rotational velocities of the swimmer for different values of C 0 and C 1 and a sideways bead attachment of x b = 0 and y b = −r, as shown schematically in Fig 7B. The full analytic forms of α i and a 0 i (i = 1, 2, 3) are very long and not particularly informative, so we only present closed form expressions for the coefficients α i , defined by Eqs 32-34, when C 0 = 0 in the supporting information, see Section 3. In particular, Eq 22 shows the closed form of α 1 . For α 2 and α 3 , Eqs. S23-S24) show the expressions for the even more restricted case where C 0 = 0 and y b = −r.
We also performed numerical simulations to study the effect of an asymmetric bead-axoneme attachment on the swimming motion (see Fig 7). In these simulations, the model  Video. An asymmetric bead-axoneme attachment causes the axoneme to rotate, as illustrated by Fig 7B and 7C; see also the S4 and S5 Videos. Thus, consistent with Eqs 32-34, our numerical simulations show that an asymmetric attachment of the axoneme to the bead causes rotation of the swimmer, even in the absence of static curvature (C 0 = 0).
It is interesting to compare the effect of the intrinsic curvature (Eq. S16 with r = 0) versus asymmetric bead attachment (Eq 32 with C 0 = 0) on the rotational velocity of the swimmer. To this end, Fig 7D presents   hO z i is comparable to that of the intrinsic curvature, C 0 . We also observe that, as shown in Fig  7E, the maximum rotational velocity is found at values of y b close to −r.

Analysis with the experimental waveform
To confirm that the predictions of the previous subsection of the existence of an anomalous propulsion regime and of a rotation induced by asymmetric cargo attachment is general and not limited to the very simplified waveform given by Eq 19, we also used the experimental beat patterns and performed RFT simulations to compute mean translational and rotational velocities of an axonemally-propelled bead for both asymmetric and symmetric bead-axoneme attachment and various bead radii.
For this purpose, we used the experimental beat pattern shown in Fig 3A of Ref. [14] (see S6 Video). As explained in our recent studies [14,22], we performed principal component analysis (PCA) of the experimental beat patterns, and then, decomposed the eigenmodes as Fourier series. Our analysis reveals that the traveling curvature waves can be decomposed into a static component C 0 and a leading traveling wave component of amplitude C 1 that coexist with standing waves at the traveling wave number and at multiples of this wave number (higher harmonics). This Fourier analysis of the experimental data indeed justifies the decomposition of the waveform as defined in Eq 19 for our analytical study.
The results of our determination of the mean translational and rotational velocities as a function of the bead radius is shown in Fig 8. Although the results are quantitatively different from those in Fig 5 obtained with the simplified waveform given by Eq 19, the variations of hO z i and one of the translational velocity show a non-monotonic dependence on r = R/L for the case when the attachment is symmetric, as shown in Fig 8A. The results obtained with the experimentally realistic waveform are therefore qualitatively consistent with those obtained with the simplified waveform, Eq 19.
In the case of an asymmetric attachment, depending on the sign of the static curvature of the axoneme C 0 and the position at which the bead is attached, one may observe an increase or decrease in the overall mean rotational velocities. For the axoneme in Fig 8, C 0 is negative and the sideways bead attachment at Y B = R (Fig 8B) acts against the rotation induced by the intrinsic curvature C 0 . The opposite happens in Fig 8C where the bead attachment at Y B = −R amplifies the rotational velocity of the axoneme. We also note that the anomalous propulsion regime is more pronounced in panel C where the bead is attached sideways at Y B = −R.
Furthermore, in the right panels of Fig 8A-8C, the measured values of hO z i/ω 0 using RFT at the experimental bead size of R = 0.5 μm (so the ratio R/L = 0.05), are indicated by red circles. We note that in this experiment (see S6 Video) the axoneme globally rotates around 2π in the time interval of *650 msec, and with f 0 = 38.21 Hz, results to hO z i/ω 0 *0.04, which is larger than the measured RFT values of 0.025, 0.020 and 0.029, corresponding to different bead-axoneme attachment geometries in panels A-C, respectively. Overall, we observe a semiquantitative agreement between the RFT predicted and the experimentally measured values of hO z i/ω 0 .
An important conclusion in Ref. [14] is that at higher calcium concentration, the mean curvature C 0 is strongly reduced, leading to a strong reduction of hO z i. For this reason, we also considered the beat patterns from the experiment in Fig 3D of Ref. [14] at higher calcium concentration (S7 Video), to study the influence of the flagellar waveform on the propulsion of the swimmer, both with symmetric and asymmetric bead attachments. The results, presented in Fig 9 also demonstrate the existence of an anomalous propulsion regimes as highlighted by the bands in cyan color in the three graphs on the right. The experimental value of hO z i/ω 0 * 0.004 (total rotation of * π/4 in 1299 msec; see S7 Video) is slightly larger than the values 0.0024, 0.0029 and 0.0028 in panels A-C, respectively, which are highlighted by the red circles in Fig 9A-9C. Finally, comparing Figs 8 and 9 shows that, as expected, the dependence of the mean translational and rotational velocities on the bead size is highly sensitive to the flagellar waveform.

Conclusions
In this work, we have studied analytically and by numerical simulations the motion of a bead propelled by a model flagellum. We used data from our previous experimental study in which isolated and demembranated flagella of the green alga C. reinhardtii were reactivated with ATP to propel a bead. [14]. In this work, we observed two distinct regimes of bead propulsion depending on the calcium concentration. The first regime describes the bead motion along a curved trajectory which is observed in experiments at zero or very small concentration of calcium ions (less than 0.02 mM). In the second regime and at higher calcium concentrations, the cargo is propelled along a straight trajectory, at an averaged velocity as high as *20 μm/sec, comparable to the typical human sperm migration speed in mucus [36]. Calcium ions are known to affect the flagellar waveform by reducing the mean curvature (C 0 ) of axonemes in a dose-dependent manner [22,23], thereby inducing a transition from circular to straight swimming trajectories.
To characterize the motion, we first used a simplified waveform to describe the axonemal shapes which is composed of a traveling wave component propagating along a circular arc (Eq 19). This simplified waveform allows us to obtain analytical expressions for the translational and rotational velocities of an axonemally-propelled bead in the limit of small amplitudes of curvature waves. The rotational velocity of an axoneme is predominately controlled by its mean curvature C 0 . As shown in Ref. [19], the second harmonics (as well as higher harmonics of even order) of the flagellar waveform also contribute to the rotational velocity of an axoneme, although more weakly (at higher orders). Remarkably, our analysis with the simplified waveform predicts a non-monotonous dependence of the rotational velocity, and/or of some of the components of the translational velocity as a function of the size of the bead. Namely, some of these components may increase when the size of the bead, hence the overall drag, increases, see Fig 5. It is also very interesting to note that the translational velocity components U x and U y are nearly saturated for a fairly large range of cargo size.
Further, we used our experimental beat patterns from Ref. [14] to demonstrate that this counter-intuitive regime is not limited to the simplified waveform and also exists for waveforms closer to the experimental ones. This anomalous propulsion regime has also been predicted for a model sperm-like swimmer with a zero mean curvature, propelled by a traveling wave component. Consistent with this, we also observed anomalous regimes using our experimental beat patterns from Ref. [14] at an increased calcium concentration (0.1 mM instead of 0 mM), in which the mean curvature of axonemes is significantly reduced (by a factor of about 10). An anomalous cargo transport regime was also predicted in biofilm forming bacteria Pseudomonas aeruginosa (PA14) [12,13], where swimming is driven by multiple (on average two) rotating helical flagella (length*4 μm) that can bundle to propel the bacterium in a corkscrewlike motion or unbundle to change direction, exhibiting a run-and-tumble swimming pattern [37]. This anomalous behavior is expected to exist for a hypothetical mutant of PA14 which has a larger (around three times) size than the wild type. This up-scaling of bacteria size results in a larger rotational drag coefficient of the flagellum, compared to that of the bacterial body, which in Ref. [12] appears to be as criterion for anomalous propulsion. Whether such a hypothetical large-scale bacterium exists is an open question. However, this anomalous propulsion regime could be important in bacterial swimming in polymeric solutions, where due to steric interactions between flagella and polymers, the rotational drag coefficient of the flagella can become larger than that of the bacterial body, fulfilling the criteria for anomalous propulsion. Thus, as experimentally observed and contrary to our expectations, the swimming speed of bacteria in the polymeric solutions can increase [38,39]. In our system, however, we are not able to obtain a simple analytical criteria for the anomalous propulsion as it results from the full calculations which include inverting the full time-dependent, 3 × 3 drag matrix, D A + D B , see Eqs 16 and 17 of the bead-axoneme swimmer. From a general physics perspective, we remark that the anomalous regime corresponds to a change in the partitioning between translation (in the two physical directions) and rotation as R increases, so that some components may increase with R over a range while other components decrease or remain almost constant over a fairly large range of the cargo size, as imposed by the overall increase of the drag due to the bead.
Furthermore, our analysis shows that asymmetric cargo-axoneme attachment provides a contribution to the rotational velocity, comparable to that of the mean curvature of the flagellum. In other words, a sperm-like beating flagellum without mean curvature and second harmonic swims in a curved trajectory if it is attached sideways to a cargo. In our experiments [14], the limitations due to the 2D imaging technique prevented us from precisely distinguishing symmetric versus asymmetric bead-axoneme attachments. Indeed, in a 2D-projected image, a symmetric bead-axoneme attachment could in reality be an asymmetric one. Moreover, as the bead-axoneme swimmer goes slightly out of focus, the attachment in some frames seems to be symmetric and in other frames asymmetric. The 3D microscopy techniques utilized in Ref. [25] are necessary to distinguish a symmetric from an asymmetric axoneme-bead attachment. This 3D characterization is absolutely essential to experimentally prove the anomalous behavior predicted by our analysis with a simplified waveform as well as with the experimental beat patterns. Although in Ref. [14], we performed few experiments with beads of diameters of 1, 2 and 3 μm, we are unable to verify the validity of the predicted anomalous trend because the 2D microscopy does not allow us to distinguish between symmetric and asymmetric bead-axoneme attachment.
Finally, it is important to note that in our analysis we have assumed that the presence of the bead (load) does not affect the waveform of the flagellum. S4 Fig in S1 File shows the hydrodynamic force distribution along the contour length of a flagellum with a given waveform at a fixed time and for different bead radii, and indeed we see that the force distribution depends on the size of the bead. These tangential and perpendicular components of the hydrodynamic force (which vary with the bead size) might feedback into the activity of the dynein molecular motors and thereby, change the flagellar waveform. By taking the product of the force distribution with the velocity distribution and integrating, we find the hydrodynamic energy dissipation to range between 0.01 × * 10 −15 J/s (r = 0) and 0.17 × 10 −15 J/s (r = 0.5). From the energy budget point of view, the ATP consumption measurements at the single-axoneme level [15] show that the energy required to generate elastic deformation in an axoneme is one order of magnitude larger than hydrodynamic dissipation. The WT active Chlamydomonas axonemes consume approximately 10 6 ATP molecules/s, corresponding to an energy consumption of 8.1 × 10 −14 J/s [15]. It is estimated that this energy is expended primarily for elastic deformation of the flagellum and not for overcoming viscous drag, as they calculate the viscous losses as 6.4 × 10 −15 J/s [15]. This suggests that the feedback effect of hydrodynamic drag on motor activity might be negligible. In-depth studies to incorporate the feedback between motor activity and hydrodynamic forces require a microscopic description of flagellar dynamics and are the subject of our future work.
The design and fabrication of synthetic micro-swimmers is a challenging task in the growing field of smart drug delivery, and has recently become a multidisciplinary effort involving physicists, biologists, chemists and materials scientists. Our theoretical analysis as well as numerical simulations reveal the existence of an anomalous cargo transport regime, where contrary to expectation, the flagellar-propelled cargo rotates faster as we increase the cargo size. This counter-intuitive behavior may play a crucial role in the design of future artificial flagellar-based propulsion systems, where targeted transport of cargo is the goal and higher rotational speeds could reduce the efficiency of directional propulsion. Finally, our analysis also highlights the contribution of the asymmetric cargo-flagellum attachment in the rotational velocity of the micro-swimmer. This turning mechanism should be also taken into account in manufacturing bio-inspired synthetic swimmers where a directional targeted motion is critical for delivery of drug-loaded cargoes.
Supporting information S1 File. Supplementary material to the manuscript. (PDF) S1 Video. Experiment: An axonemally-driven bead. An isolated and demembranated flagellum from green algae C. reinhardtii is attached to a 1 micron-sized bead. The axoneme is reactivated with 1 mM ATP and beats at around 110 Hz.